Variable VHE 7-ray emission from Markarian 501 
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ABSTRACT 

The blazar Markarian 501 (Mrk 501) was observed at energies above 0.10 TeV 
with the MAGIC telescope from May through July 2005. The high sensitivity of 
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the instrument enabled the determination of the flux and spectrum of the source 
on a night-by-night basis. Throughout our observational campaign, the flux from 
Mrk 501 was found to vary by an order of magnitude. Intra-night flux variability 
with flux-doubling times down to 2 minutes was observed during the two most ac- 
tive nights, namely June 30 and July 9. These are the fastest flux variations ever 
observed in Mrk 501. The ~20-minute long flare of July 9 showed an indication 
of a 4±1 min time delay between the peaks of F(<0.25 TeV) and F(>1.2 TeV), 
which may indicate a progressive acceleration of electrons in the emitting plasma 
blob. The flux variability was quantified for several energy ranges, and found to 
increase with the energy of the 7-ray photons. The spectra hardened significantly 
with increasing flux, and during the two most active nights, a spectral peak was 
clearly detected at 0.43 ± 0.06 TeV and 0.25 ± 0.07 TeV, respectively for June 
30 and July 9. There is no evidence of such spectral feature for the other nights 
at energies down to 0.10 TeV, thus suggesting that the spectral peak is correlated 
with the source luminosity. These observed characteristics could be accommo- 
dated in a Synchrotron-Self-Compton (SSC) framework in which the increase in 
7-ray flux is produced by a freshly injected (high energy) electron population. 

Subject headings: Markarian 501, BL Lac, AGN, VHE 7-ray, imaging air Cherenkov 
telescope, MAGIC 



Introduction 



The large inferred luminosities of Active Galactic Nuclei (AGNs) led to a standard model 
of beamed AGN emission, with the ultimate energy source being the release of gravitational 
potential en ergy of matter from an accretion disk surrounding a super-massive black hole 
(Rees 1984 ). Particularly interesting for the Very- High- Energy 7-ray (VHeI^ community 
are the blazars, whose relativistic plasma jets point at the observer. Distinctive features of 
blazars are their continuum emission, clearly non-thermal from radio to VHE frequencies and 
characterized by two broad bumps peaking at, re s pectively, IR/X-ra y and 7-ray frequencies 
(IBlandford fc Reed Il978l : lUrry fc Padovani I Il995l : lUlrich et al 1 119971 ) , and their strong vari- 
ability, implying flux variations by a factor of <^ 10 over timescales of ^ 1 hour to months 
(flares; see Ulrich et al. 1997). 



So far, 13 AGNs have been detected at VHE energies. Except for M87 (lAharonian et al 



^ In this paper the VHE band is defined as the energy range E>0.1 TeV. 
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20031 . l2006l ). all these sources belong to the 'high-peaked BL Lac' (HBL) sub-class of blazars, 
which are characterized by a spectral energy density (SED) in which both maxima occur 
at relatively high frequency (e.g., respectively at hard X-rays and HE/VHE 7-rays). The 
detection of 7-rays from blazars leads to some important considerations about the relevant 
radiation processes and the physical properties of the emitting regions. Most notably, the 
very detection of VHE radiation, implying 7-ray transparency in the emitting region, re- 
quires the presence of relativistic beaming to decrease the intrinsic energy density of the 
soft target photons inside the source. The beaming reduces this value because it simultane- 
ously decreases the intrinsic energy density of the photons and, having reduced the energy 
of the relevant 7-ray photons, it increases the energy of the soft target photons relevant to 
e^ production, hence (for typical spectra) decreasi ng their number den sity (McBreen 1979; 



f lBlandford fc Konigllll979[ ). Bassani & Dean 1981; flMattox et allll993f ): Dondi & Ghisellini 
1995). 

Two classes of emission models have been proposed to explain the TeV emission from 
blazars: leptonic and hadronic models. 

(i) In the case of the most popular leptonic models, the same population of non-thermal 
electrons (and possibly positrons) responsible for the radio-to-X-rays SED is also responsi- 
ble for 7-ray emission, through Compton up-scattering of the synchro tron photons off their 



own parent electrons - the Synchrotr on-Self-Compton (SSC) process ( iMarscher et allll985 



Maraschi et al 1119921 : iBottcher II2002I ). In other models, electrons scatter 'external' photons 
that originate outside the jet (External Compton (EC) model^. In BL Lacs the lack of 
strong emission lines suggests a minor role of ambient photons and hence supports the SSC 
models. 

(a) In hadronic models, the TeV radiation is produced by hadronic interactions of the highly 
relativistic baryonic outflow with the ambient mediunj^, and/or by i nteractions of ultra-high 



energy protons with synchrotron photons produced by e lectrons (IMannheim fc Biermann 
1992!). with the jet magnetic field (Aharonian et al 2000), or with both ( Miicke et al 
Atovan et al II2OO3I : iMannheim Ill993h . 



2003 



However, hadronic models are challenged by the blazars' observed X-ray versus VHE corre- 
lation and very rapid 7-ray variability. The SSC model is then widely believed to explain 



The external seed photons may, e.g., come from the accretion disk ( Dermer fc SchUckeiser Il993h . 



be the disk radiation scattered by the material around the disk and and the jet ( [ Sikora et all 11994 ) 



be radiation from the massive stars which enter the jet (jBednarek fc P rotheore Ill997a[). or be syn chrotron 



or 
or 



radiation produced in the jet and reflected by the surrounding material (.Ghisellini fc Madau 19961 ) 



I.e.: gas and clouds drif ting across the 



et (jPar fc Laor Ill997t iBeall fc Bednarek 1999). the matter 



from the thick accretion disk ( Bednarek Ill993| ). or interactions inside the (dense) jet ( Pohl fc Schlieckeiser 
200(h . 
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the dominant emission process in blaza r jets (not always in its s imple st one-zone realizati on 



as required by e.g. 'orphan flares': see iKrawczynski et all (120041 ) and iGliozzi et al I (120061 )). 



In this framework, the importance of high-quality data on blazar VHE emission can 
not be overestimated. In particular, valuable information can be obtained investigating: (i) 
the rapid, possibly energy-dependent, flux variability; (ii) the X-ray/optical versus VHE 
correlation; and (Hi) the X-ray and VHE spectral variability, with a potential energy shift of 
the Synchrotron and Compton peaks. Simultaneous multi-wavelength observations of such 
rapid variability can provide stringent tests to emission models, in particular on acceleration 
processes in jets. In the case of nearby {z<0.1) sources, the extinction due to pair production 
by interaction of the blazar-emitted TeV photons with the (not well known) optical/IR 
Extragalactic Background Light (EBL) photons is probably minor, and thus there is a smaller 
uncertainty in the determination of the intrinsic spectral features of the object. 

In this paper we report about detailed measurements of the VHE emission of Mrk 501, 
demonstrating the capability of the MAGIC Telescope in the precision stu dy of blazar 



physics. The BL Lac objec t Mrk 501 was the second established TeV-blazar (IQuinn et al 



19961 : iBradbury et al 1119971 ). After a phase of moderate emission for about a year following 
its discovery as a TeV source, in 1997 Mrk 501 went into a state of surprisingly high activity 
and str ong variability, becomiii g >10 times brighter (at energies >1 TeV) than the Crab 
Nebula (lAharonian et al lll999al Jbl'). In 1998-1999 the mean flux dropped by an order of mag- 



nitude (lAharonian et al 



200ll ). It is worth noticing that the HEGRA observations (with 



a threshold energy of ~0.5 TeV) did not see spectral variations during the 1997 outburst, 
whereas it did observe in 1998-1999 a significantly softer low-state energy spectrum than in 
1997. The CAT telescope (with a thr eshold energy of 0.25 TeV), on the other hand, did 
detect spectral variations during 1997 ( Djanati-Atai'et al 1999; Piron 200ol . [2OO1I ) . 



The structure of the paper is as follows. In sect. 2 we briefiy describe the MAGIC data 
taking and analysis. In sects. 3 and 4 we present and discuss the Mrk 501 VHE lightcurve 
(hereafter LC), spectrum and their variability, during the observation campaign. In sect. 5 
we discuss the long-term light curve of Mrk 501, short-term fiux variability, fiux-spectrum 
correlations, and the overall SED of this object. Finally, sect. 6 summarizes our main results. 



2. The MAGIC Telescope and the data analysis 



2.1. The instruments 



The observations in the VHE domain were carried out with the Major Atmospheric 
Gamma-ray Imaging Cherenkov (MAGIC) telescope, located on the Canary island of La 
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Palma (28.8''A^, 17.9°!^) at the Roque de los Muchachos Observatory (about 2200 m above 
see level). MAGIC started regular observations in the fall of 2004 and, with a main mirror 
diameter of 17 m, it is currently the world's largest single-dish lACT. Further details about 



the characteristics and performance of MAGIC c an be found elsewhere ( iBaixeras et al 112004 
Paneque l[20oi : ICortina et aTlboOsI : ICaug Ibooeh . 



The MAGIC Collaboration also operates the optical KVA telescope (35 cm). Simulta- 
neously with the MAGIC observations, Mrk 501 was regularly observed with KVA as a part 
of the Tuorla Observatory blazar monitoring progranQ. In this paper we also use 2-10 keV 
data taken with the RXTE satellite's All-Sky-Monitor (i?XrE/ASM)i quasi- simult aneously 
with our MAGICobservations. 



2.2. Source observation 

The source was observed during 30 nights between May and July 2005, with an overall 
observation time of 54.8 hours. In order to maximize the time coverage of this source, 
observations were carried out also in the presence of moonlight (34.1 hours, i.e. 62% of 
the total observing time). It is important to note that many of these 'moon observations' 
were performed when the moon was only partly illuminated, and mostly located at a large 
angle (60°-90°) with respect to the position of Mrk 501. This kept the Night Sky Background 
(NSB) of these observations rather low and comparable to that of the moonless observations. 
The observations were mostly performed in the so-called 'on mode' in which the telescope 
points exactly to the source {on- data), and thus its (optical) image is right in the center of 
the camera. The telescope was also operated in 'off mode', in which it points to regions of the 
sky where there are no known 7-ray sources {off-data). These observations were carried out 
at comparable zenith angle (ZA) and NSB conditions, and can therefore be used to estimate 
the background content of the on-data. The off-data observation time amounts to 3.5 hours. 

The data were screened for hardware problems, non-optimal weather conditions, and 
too high NSB light. In addition, the very few runs with ZA larger than 30° that survived 
those filters were also removed, in order to have a more uniform data set. The number of 
nights surviving the=e .election cut= i= 24, with a total net observation time ot 31.6 howfl 



http: //users .utu. f i/kemi/lm/ 



The data are publicly available at http://heasarc.gsfc.nasa.gov/xte_weather/ 



^ Five out of the six entirely rejected observation nights (with a total observation time of 9.6 hours), were 
moonlight observations. 
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These observation nights, together with the corresponding net times and ZA ranges of the 
data acquisition, are hsted in Table [H 



2.3. Data analysis 



The analysis used in this paper is based on the Hillas image parameters A LPHA, 
WIDTH, LENGTH, DIST and SIZE to quantify the camera image, see lmUasI JlQSsh . These 
parameters are calculated using the calibrated signals from the individual pixels of the cam- 
era. T h e pro cedur es used to cal i brate the PMT signals in MAGIC are described in detail in 
Gang I (120061 ) and lAlbert et al I (120061 ): here we used the Sliding Window signal extractor, 
and performed the calibration with the Excess Noise Factor method. Only PMT signals 
with more than 10 photoelectrons (8 photoelectrons for pixels in the boundary of an image) 
which occur within a time window of 6.6 ns (2 FADC slices) of the neighbouring pixel signal 
were used. The minimum total image light content (SIZE) considered in this analysis is 150 
photelectrons. The signal/background separation is achieved by applying dynamical cuts (de- 
fined as 2nd order polynomial functions of logSIZE) on the parameters WIDTH, LENGTH 
(shape parameters) and DIST (position of the image). The background contained in the 
on- data after the 7/hadron separation cuts is estimated by means of a 2nd-order polynomial 
fit (with no linear term) to the ALPHA distribution from the normalized (according to the 
on-off observation times) off-data. 

The energy of the incoming 7-rays is essentially proportional to the light content of the 
image (SIZE), with corrections according to the values of LENGTH, DIST and LEAKAGE^. 
The energy resolution achieved with this parameterization is about 20-30%, slightly depend- 
ing on the event's energy. Because of the finite experimental resolution, the distribution 
of the excess events (7-candidates) versus the reconstructed energy is a convolution of the 
(true) energy distribution of excess events and a realistic energy resolution function. The 
determination of the true e nergy spectrum from the reconstructed one is achieved by means 
of an unfol d ing p rocedure (lAnykeyev et allll99ll ). We used the iterative method described 



m 



Bertero I (119891 ). 



In this analysis, 90-99% of the background images are removed by the selection cuts, 
while 50-60% of the 7-ray signals are kept. The resulting collection area after analysis cuts 
is 0.5 • 10^ down to 0.20 TeV. The analysis threshold energy, commonly defined as 
the peak of the differential event-rate spectrum after all cuts, is ~ 0.15 TeV. The lowest 



^ LEAKAGE is defined as the fraction of the hght content recorded by the outer ring of the PMT camera, 
and it is typicaUy used to evaluate the level of missing light in the detected image. 
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7-ray energy used in the calculation of the energy spectra is 0.10 TeV. In the LC analysis, 
however, the minimum energy considered is 0.15 TeV (i.e. the threshold energy). Below 
this energy the collection area drops fast, increasing rapidly both systematic and statistical 
errors in the measured flux. Keeping the measurement errors small is essential for the study 
of flux variations, one of the main goals of this work. 

In order to check the reliability of the used analysis chain, we analyzed data from 
the Crab Nebula taken in December 2005, under similar instrumental and environmental 
conditions to those of Mrk 501. The obtained results were in perfect agreement w ith data 
published previously (IHillas et al Ill998l : lAharonian et al Il2004j : IWagner et al 1 120051 ) . which 
shows that the analysis procedures used produce reliable results. 

The results from the MAGIC Mrk 501 data analysis are summarized in Tabled! The 
table shows the integrated flux (above 0.15 TeV) and the resulting fit to the differential 
(energy) photon spectra with a simple power-law (PL) model (see eq. |5]), for each observing 
night. The fit was obtained using all the spectral points above 0.10 TeV. Only statistical 
errors are reported in Table [U The systematic errors on the energy determination are esti- 
mated as ~20% which, for a spectral index of 2.5, would produce a systematic shift of 50% in 
the flux level (normalization factor of the PL function from Tabled]). The systematic error 
in the calculated spectral indices is evaluated as ~0.1. In the tab le we quote the c ombin ed 
significance which is calculated following the prescription given by iBityukov et al I (120061 ) as 



S, 



comb 



where Si is the significance corresponding to the (differential) energy bin i, and n is the 
number of energy bin s (measurements ). The significance of each energy bin is calculated 
according to eq.(5) of iLi &: Ma I (119831 ) . which is more suitable than eq.(17) from the same 
paper, since in all nights the 7-ray signal is clear and hence its existence is not in doubt. 
The combined significance is used to compare the quality of the 7-ray signals from different 
observing nights. 



3. Lightcurve of Mrk 501 during the MAGIC observations 

In this section we report on the broadband, optical to 7-ray, LC of Mrk 501 during 
May- July 2005. 
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3.1. Overall lightcurve at 7-ray, X-ray, and optical frequencies 

The overall LC of Mrk 501 during the MAGIC observation campaign is shown in Fig. 
[U The observed flux is shown in three energy bands: VHE (0.15TeV-10TeV), X-rays (2keV- 
lOkeV), and optical (1.5eV-2.5eV) as measured by MAGIC, RXTE/ ASM and KVA, respec- 
tively The X-ray and optical fluxes are computed as weighted averages using RXTE/ ASM 
and KVA measurements taken simultaneously with the MAGIC observations plus/minus a 
time tolerance of 0.2 days. A smaller time tolerance substantially decreases the number of 
X-ray points that can be used. The flux level of the Crab Nebula (lilac-dashed horizontal line 
in the top plot) is also shown in Fig. [T]for comparison. The Crab Nebula flux was obtained 
by applying the very same analysis as described in Sect. 12.31 to the MAGIC Crab Nebula 
data taken during December 2005 under observing conditions similar to those for Mrk 501. 
The estimated Crab Nebula flux level is therefore roughly affected by the same systematics 
as the fluxes obtained for Mrk 501. We found Fcrab(>0.15 TeV) = (3.2±0.1) x 10"^° cm'^s'^, 
thereafter referred to as Crab Unit (c.u.). For simplicity, only the Crab Nebula flux level, 
and not the associated error (which is irrelevant for the comparison) is shown in the LC. 

The measured VHE flux from Mrk 501 was at about 0.5 c.u. during most of the observa- 
tion nights (Table [1]). During several nights, however, its flux signiflcantly exceeded 0.5 c.u., 
and during one night (MJD 53536.947) it showed a substantially lower flux (0.24±0.04 c.u.). 
Often, Mrk501 showed large flux variations in consecutive nights. An example of these rapid 
flux variations are the MJD 53535.934 and 53536.947 with respective fluxes of 0.84±0.04 c.u. 
and 0.24±0.04 c.u.; the MJD 53554.906 and 53555.914 with respective fluxes of 1.11±0.09 c.u. 
and0.40±0.11 c.u.; and the MJD 53563.921 and 53564.917 with respective fluxes of 1.74±0.09 c.u. 
and 0.91±0.15 c.u.. Besides, the VHE flux from Mrk 501 was outstanding during the MJD 
53551.905 (June 30) and 53560.906 (July 09) with 3.48±0.10 c.u. and 3.12±0.12 c.u., respec- 
tively. During these two nights the source was in a very active state. Note, however, that 
the night before the July 9 flare the emitted flux was 0.58±0.07 c.u., i.e. close to the average 
flux of the entire campaign. Mrk 501 therefore showed a remarkably fast VHE variability 
during this campaign. 

Unlike in VHE 7-rays, no signiflcant flux variation was recorded in the X-ray and 
optical bands. In the case of the X-ray data, however, although the sensitivity of the 
i^XTi?/ ASMinstrument was clearly inadequate to reveal short-term 2-10 keV flux variability 
in Mrk 501's emission, the flux appears to be higher in the second portion of the LC. The 
optical flux, on the other hand, shows only a modest variation, a ~5% monotonic decrease 
during the entire observational campaign. 
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3.2. Multi-frequency correlations 

The correlations of our observed VHE 7-ray data with X-ray and optical data are shown 
in Fig. [21 the gamma points being the same as shown in the light curve of Fig. [TJ It can be 
seen that the measurement uncertainties of the X-ray and optical fluxes are comparatively 
large, which makes possible different dependencies of the X-ray flux on the 7-ray flux. We 
flt the VHE/X-ray data with a constant and with a linear function, obtaining the highest 
probability for the linear function. For the relation VHE/optical data, the two flts are 
nearly equally probable. All flt results are shown in the insets of Fig. [21 including for the 
VHE/X-ray data also a flt forcing the linear function through the origin. 

For the VHE/X-ray correlation, one obtains a linear correlation coeflicient of 0.49. In- 



( 


2007 


); 


Hrupec 


( 


2007) 



points (corresponding VHE and X-ray fluxes), multiple (we used 100 000) possible sets of 
measurements are generated, using random differences derived from the (Gaussian) errors 
assigned to each point. For each set of generated points, a correlation coefficient is obtained, 
resulting for a large number of measurements in the probability density function (pdf) of 
correlation coefficients /a('") shown in Fig. [3l which corresponds to the measured correlation 
and the assigned errors. This pdf correctly expresses the effect of the measurement uncer- 
tainties. The same procedure can be applied to hypothetical fully correlated data: to this 
end, the data points are shifted onto the straight line from the fit with the highest prob- 
ability (black line in left-hand plot of Fig. [3]), maintaining the original error assignments. 
In the case of hypothetical uncorrelated data, the points are randomly distributed in the 
VHE/X-ray plane. Sets of Monte Carlo measurements are then generated as before. The re- 
sulting pdfs express the probability to obtain certain values of correlation coefficients, given 
respectively no correlation (/c(r) in Fig. [3]) or full correlation (/^(r) in Fig. [3|) with our 
error assumptions. For comparison purposes. Fig. [3] also shows the analyti cal pdf functio n 



for the uncorrelated case, foif) (described in chapter 9 and appendix C of (iTaylor Ill997l )). 
which does not take into account measurement errors. Note that /c('") is very similar to 
foif)- This is not surprising, since the smearing does not increase the randomness of the 
already randomized seed event. The width of the probability density distribution foir) de- 
pends exclusively on the number of points in the data sample (23 in our case). On the other 
hand, /^(r) is significantly affected by the measurement errors; even under the hypothesis of 
a fully correlated case, the probability to obtain values for the correlation coefficient larger 
than 0.8 is very small. 

A measure of the probability of correlation can be derived from the comparison of the 
probability density distribution /a('") for the actual measurement with the distributions for 
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the two extreme correlation cases, /^(r) and fcij')- Given the results in Fig. [3l it is evident 
that is similar to /^(r), but, despite a sizable overlap, rath er different frorn fci^)- 



For a quantitative comparison, we followed the robust method from lPoe et all (120051 ) . which 
is based on the convolution of empirical probability density distributions. The comparison 
of a pair of probability density functions fx{r) and /y(r) leads to the probability that 
these two distributions are statistically consistent, P{fx{r),fY{r)). The resulting value 
for the probability of agreement between our data points and the fully correlated case, 
Px,'y{fA, fs), is 0.55 ± 0.15, being the quoted error a systematical error estimated through 



variations in t he initial values used in the Monte Carlo method (IFerenc fc Hrupec 1 12007 



Hrupec 1 120071 ). The probability that our measurement came as a result of a statistical 
fluctuation of an entirely uncorrelated physics case, on the other hand, is significantly lower; 
Px,'r{fA, fc) is 0.15±0.05. The probability for the first scenario to be true and for the second 
to be false is Px,-y{fA,fB){^ — Px,'y{fA, fc)) = 0.47, and the probability for the opposite 
is (1 - Px, ■y{fA,fB))Px^{fA,fc) — 0.067, which indicates that the correlation scenario is 
significantly more likel\f^. 

The same method was applied to the optical and 7-ray flux values of Fig. [2] (right-hand 
plot). The corresponding probability density functions are presented in Fig. |H In this 
case, the resulting linear correlation coefficient is -0.27, indicating a small anti-correlation; 
yet the probabilities -Popt,7(/yi, /b)=0.60±0.25, and -Popt,7(/A, /c)=0.55±0.25 are practically 
equal, which suggest that neither the anti-correlated, nor the uncorrelated scenario may be 
reliably excluded. The probability for the first scenario to be true and for the second to be 
false is 0.27, while the probability for the opposite case is 0.22, which confirms the previous 
conclusion. 



3.3. Intra-day 7-ray flux variations 

During the two nights with the highest VHE activity, namely June 30 and July 9, Mrk 
501 clearly showed intra-night flux variations. The corresponding LC in the 0.15-10 TeV 
band is shown in Fig. [5] with a time binning of ~2 minutes. For comparison, the Crab 
Nebula flux is shown as a lilac dashed horizontal line. The vertical dot-dashed line divides 
the data into a region of relatively 'stable' (pre-burst) emission and one of 'variable' (in- 
burst) emission. The background rate after the gamma/hadron selection cuts was evaluated 
during these two nights, and is shown in the bottom panels of Fig. O These rates were 



^ It is worth noticing that the present analysis treats data from low and very-high activity epochs together, 
whereas the same correlation slope may not be necessarily the same for the two activity states. 
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found to be constant along the entire night. Consequently, the variations seen in the upper 
panels of Fig. [5] correspond to actual variations of the VHE 7-ray flux from MrkSOl, thus 
ruling out detector instabilities or atmospheric changes. 

A constant line fit to the whole LC gives a x^/NDF = 47.9/30 (probability p = 2.0 x 
10-2) for the night of June 30, and a x^/NDF = 80.6/21 {p = 6.4-10-9) for the night of July 
9. The emission above 0.15 TeV during the two nights is therefore statistically inconsistent 
with being constant and it is more reasonable to only fit the first part, i.e. the 'stable' 
emission, with a constant - and not the entire LC. A constant fit to the 'stable' portion of 
the LCs gives x^/NDF = 13.4/12 (p=0.34) for June 30, and x^/NDF = 17.8/11 (p=0.09) 
for July 9 (see Fig. [5]). The probability that the 'variable' parts of the LCs are compatible 
with the 'stable' flux level is given by x^/NDF = 34.5/18 {p = 1.1 x 10-^) for June 30 
and x^/NDF = 83.3/10 {p = 1.1 ■ 10^^^) for July 9, respectively. We therefore measured 
intra-night flux variations in both nights. 

The flare's amplitude and duration, as well as its rise/fall times, can be quantified 
according to 

F{t) = a + (2) 
2- c +2^ 



(jSchweizer II2004I ). This model parameterizes a flux variation (flare) superposed on a stable 
emission: F{t) asymptotically tends to a when t ±00. The parameter a is the assumed 
constant flux at the time of the flare (cf. the horizontal black dashed lines in Fig. [5]); is set 
to the time corresponding to the point with the highest value in the LC; and b, c, d are left free 
to vary. The latter two parameters denote the flux-doubling rise and fall times, respectively, 
and can be converted into the characteristic rise/fall timecl by multiplying them by l/ln2. 
The resulting fits using eq. [2] are shown in Fig. [5|, and their parameter values are reported 
in Table [21 In both cases, the measured rise/fall flux-doubling times are ~2 minutes, which 
yields characteristic rise/fall times of 2/ln2 ~3 minutes. These are the shortest flux- variation 
timescales ever measured from Mrk 501, at any wavelength. 

Because of the steeply falling spectra, the low-energy events dominate the LCs shown 
in Fig. [5l and tend to hide any higher-energy features. We therefore split the data into 
four distinct energy ranges: 0.15-0.25 TeV, 0.25-0.6 TeV, 0.6-1.2 TeV and 1.2-10 TeV. The 
corresponding LCs for the night June 30 are shown in Fig. O Due to the reduced photon 
statistics, we increased the time binning from 2 minutes to 4 minutes. We found that 
only the energy range 0.25-0.6 TeV shows a clear flux variation; a constant line fit gave 
X^/NDF = 24.2/8 [p = 2 ■ 10-^). The other energy ranges are compatible with a constant 



The characteristic time is defined as the time needed for the flux to change by e^^. 
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line fit, showing only a slight overall flux level variation with respect to the LC in the 'stable' 
part. Therefore, if there is a flux variation, it is too small to be signiflcantly seen in our data. 
Note that, among the 4 energy ranges used, the 0.25-0.60 GeV energy range is the one with 
the highest sensitivity for flux variations. The 'variable' LC from the energy range 0.25-0.60 
TeV was then fltted using equation [2], flxing to to the time of the highest point in the LC. 
The resulting parameters of the flt are reported in tabled the rise/fall flux-doubling times 
are comparable to those ones obtained using the integrated LC above 0.15 TeV. 

The same exercise on the flare July 9 gave a signiflcantly different result, as shown in Fig. 
[71 The flare is visible essentially in all energy ranges. In order to study possible time shifts 
between the different energy ranges, we flt all LCs simultaneously (combined flt) with a flare 
model described by equation [2J In order to remove one degree of freedom and facilitate the 
flt procedure, we assume a symmetric flare with equal rise and fall flux-doubling times; that 
is c = c? in equation [21 The resulting parameters from this combined flt are shown in tabled 
The combined flt gave /NDF = 14.0/12 {p = 0.3), which implies that the measured flare 
is compatible with being symmetric. The rise/fall flux-doubling time is about 2 minutes for 
all the energy ranges. It is interesting to note that the position of the peak of the flare for the 
different LCs seems to vary somewhat with energy. The time difference between the highest 
energy range and the lowest energy range is 239 ± 78 seconds. If instead, the energy range 
0.25-0.6 TeV is selected as the lowest energy range, which has a better deflned flare (and 
thus a better determination of the peak position), the time difference is 232 ±54 seconds. 

In order to evaluate the signiflcance of this time shift, we performed the same flt, but 
this time using a common to for all LCs. The resulting flts are shown in Fig. [HI and the 
resulting parameters from the flt in table [51 The combined flt gave /NDF = 26.6/15 
{p = 0.04), which implies that such a situation is unlikely, and consequently that the time 
shift of 4 ± 1 min between the highest and the lowest energies is more probable. 

Investigating the reliability of the time delay obtained from the combined flt, we per- 
formed a cross-correlation analysis on the LCs from July 9 with the methodology described 
in section 13. 2[ For this stud y w e used LCs with 2 min time bins from the energy ranges 
0.25-0.6 TeV and 1.2-10 TeV^j. The correlation coefficient and probability of correlation 
were computed after introducing time shifts of 2 min (one bin in the LCs). We obtained 
the highest values for a time lag of 4 minutes; which is consistent with the results from the 
combined flt shown above. 



The flare observed in the LC from the energy range 0.15-0.25 TeV is not very well defined because of the 
somewhat larger measurement errors, and the smaller relative amplitude of the flux variation (with respect 
to the 'stable' emission) with decreasing energy. 
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We want to point out that this is the first time that a possible time delay between fiares 
at different energies is observed at VHE 7-ray energies, although such time lags have been 
detected for so nie TeV blazars at X -ray frequencies, viz. Mrk 421 (IRavasio et al 1 120041 ) and 



PKS 2155-304 fIZhang et al ll2006al Jbll If the observed VHE time lag is assumed real, this 
suggests that we are observing the underlying dynamics of the relativistic electrons in both 
the synchrotron and IC emission, and the observation, therefore, supports SSC models. 

It should be also noted that the relative amplitude of the fiux variations observed in 
the LC for July 9 with respect to the baseline emission, is significantly larger at the highest 
energies. This can be seen from the ratio 6/a, where a and h are the coefficients in eq. [2] 
describing, respectively, the baseline and amplitude of the fiare (see table H]): h/a = 3.6 ± 1.0 
and 17±4 for, respectively, the 0.25-0.6 TeV and 1.2-10 TeV bands. The July 9 LC also 
shows some significant fiux variation in its 'stable' part: in the highest-energy band, where 
activity is most conspicuous, a constant fit gives a /NDF = 20.6/5 {p = 9.6 ■ 10^"^). 

In summary, during the 2005 MAGIC observations of Mrk 501 we detected variability 
at VHE frequencies with fiux-doubling times down to 2 minutes. This is about 50 times 
faster than the shortest previously observed variability-times at VHE frequencies for Mrk50 1 



(iHayashida et al II1998I : iQuinn et al II1999I : lAharonian et al Ill999al : iDjanati-Atai'et al 
and about 5 times shorter than the shortest observed variability for Mrk 421 (jCaidos 



1999) 



1996|). 



The above presented fiu x- variations are among the shortest ever observed in blazars, see also 
Aharonian et al I (120071 ). It is interesting to note that the Mrk 501 fiux-doubling rise times 
observed by MAGIC in the VHE range are rather compa rable to the shorte st variability 
times observed at X-ray frequencies which were reported by IXue fc Cui I (120051 ): a fiare with 
a total duration of 15 minutes with a fiux variation of 30%. The authors, however, reported 
the presence of substructures, which point to the existence of variability on timescales shorter 
than 15 minutes. It is worth mentioning that for both X-ray and 7-ray the shortest fiux 
varitions occurred when the source was not in an exceptionally high emission state. 



3.4. Quantification of tiie variability 



Mrk 501 has shown energy- dependent fiux variations throughout the e ntire MAGIC ob- 
servational campaign. We followed the description given in lVaughan et al I (120031 ) to quantify 
the fiux variability by means of the fractional variability parameter F^ar, as a function of 
energy. In order to account for the individual fiux naeasur ement errors {(Terr,i), we used 
the 'excess variance' (INandra et al 1119971 : lEdelson et al 1120021 ) as an estimator of the intrin- 
sic source variance. This is the variance after subtracting the expected contribution from 
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measurement errors. For a given energy range, the F^ar is calculated as 



< F^ > ' 



where < > is the mean photon flux, S the standard deviation of the flux points, 
and < al^j. > the average mean square error, all determined for a given energy bin. The 
uncertainty on Fyar is estimated according to: 




Fig. [9] shows the derived F^ar values for 5 logarithmic energy bins, spanning from 
0.14 TeV to8 TeV0. The left-hand plot includes all data, while the right-hand plot includes 
all data except for the two active nights. The result indicates a larger- amplitude flux variation 
at higher energies, which is clearly discernible even when the active-state data are excluded 
(the null probability being p~10~^). 

Fig. [ini shows the F^ar values as derived for the individual active nights of June 30 
and July 9. The flux variability in these nights is smaller than for the entire observational 
campaign, as one would expect from simple inspection of Fig. [TJ In the night of June 30, 
Fyar does not increase signiflcantly with energy. In contrast to June 30, there is a clear 
increase of F^ar with energy in the night of July 9 - in spite of the larger error bars coming 
from a shorter observation time and a lower mean flux. 

In summary, during the year 2005 MAGIC observations, the VHE 7-ray flux variability 
of Mrk 501 was found to signiflcantly increase with energy, on time scales from months to 
less than an hour. A similar effect ( on tim escales ^1/2 hr) was also detected in X-rays in 



1997, 1998, and 2000 (IGhozzi et al I fl2006f ). ba sed on RXTE data) Another X-ray evidence 



was found for the TeV blazar: PKS 2155-304 (jZhang et al Il2006bl ). The largest X-ray F, 



value in X-rays for Mrk 501 was ~0.6-0.7 in the highest energy bin a t 10-20 keV, and it 



was found in June 1998 and July-September 2000 (iGliozzi et al I (l2006l )). In 1997 however, 
despite Mrk 501 showed the highest X-ray (2-10 keV) fluxes from the last ten years, the 
largest observed Fyar value in X-rays was only ~0.4. It is interesting to note that, in the 
VHE 7-ray range, we observe a maximum Fyar of ^ 0.6, which increases to ^ 1.2 when the 
two active nights are included. This suggests that Mrk 501 is more variable in VHE 7-rays 
than in X-rays. 



The Fyar is not meaningful below 0.14 TeV (i.e. below the threshold energy of the instrument) because 
the flux errors are rather large, which makes < (7^^,^ S^. 
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4. VHE Spectra 

The differential photon spectra of Mrk 501 were parameterized with a simple power-law 
(PL) function: 

dE = ^'' io^) 

where Kq is a normalization factor and a the photon index. The MAGIC sensitivity permits 
to derive VHE spectra of Mrk 501 on a daily basis, independent of its flux level, a real asset 
for unbiased precision studies of blazars. The PL spectral parameters for each single night 
are reported in Table [H 

During our observations the VHE emission of Mrk 501 showed a very dynamic be- 
haviour, with significant spectral variability on a timescale of days (see Sect. 14.11 and 14. 2p . 
Nevertheless, most of the data are well described by the simple PL function of eq. 6. This 
does not hold for the two flaring nights of June 30 and July 9 which are therefore discussed 
separately in Sect. 14.31 and I 



4.1. Spectral index vs. flux 

The Mrk 501 VHE spectrum was measured on a night-by-night basis, which allowed 
investigating possible correlations between the PL spectral index and intensity, as it is shown 
in Fig. [m The data of June 30 and July 9 were again split into the 'stable' and 'variable' 
part (see Fig. [5]) and consequently. Fig. [11] contains 26 instead of 24 points. The data points 
are well described by a linear fit; but a constant is clearly excluded {x^/NDF = 74.7/25, i.e. 
p = 7.4-10^''). On average, the Mrk 50 1 spectrum h arden s when the emission incre ases. Such 



a correlation was already reported by iPian et al I (119981 ): iTavecchio et al I (120011 ). although 
on substantially longer timescales. 

In order to test the ansatz of a linear correlation in spectral index versus 7-ray flux 
(see Fig. [TTl) . we applied a correlation analysis as described in Sect. 13.21 The linear fit 
indicates anti-correlation (x^/ndf = 19.65/24, p=0.72). The correlation factor is -0.48, and 
the corresponding probability distributions are shown in Fig. [121 We find the probability 
-Pa,7(/A, /s) = 0.35, to be higher than Pa,-y{fA,fc) = 0.12, which supports a modest anti- 
correlation. 

Similar, the probability for the first scenario to be true and for the second to be false 
is PQ,^-y(/^, /b)(1 — Pa,'r{fA,fc)) = 0.30, while the probability for the opposite case is (1 — 
Pa,'y{fA,fB))Pa,'y{fA,fc) = 0.08. In couclusiou, this Correlation study indicates a spectral 
hardening with increased flux. 
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4.2. Spectra at different flux levels 

In order to investigate the spectra at different fiux levels, the diurnal data were combined 
into three groups, depending on whether the integral flux above 0.15 TeV, -Fo.i5TeV (measured 
in Crab Units (c.u.)), was low (Fo.i5TeV<0.5 c.u.), medium (0.5 c.u.<Fo.i5TeV< 1-0 c.u.), or 
high (1.0 c.u.<Fo.i5Tev)- Based on the chosen flux limits the low-, medium- and high-flux 
data sets consist of 12, 8 and 2 nights, respectively (see tabled] for detailed statistics). The 
data from June 30 and July 9 will be discussed in Sect. 14.31 and are not included in the 
analysis here. The differential photon spectrum for all three flux regions together with the 
PL fit results are shown in Fig. [T31 the fit parameters are also listed in table O Even with 
such a simple parameterization, our data do suggest a spectral hardening with increasing 
flux. The results of this analysis are consistent with the trend seen in Fig. [11] and discussed 
in sect. 14. 1[ 

The HEGRA CT system measured the spectra of Mrk 501 in 1998-1999, at a time when 
its flux level was s ubstantially below the one of the Crab (~l/3 c.u.). The observation 



covered 122 hours (lAharonian et al Il200ll ). The reported spectrum could be fitted with a 
PL in the energy range 0.5-10 TeV, giving a spectral index of 2.76 ± O.od^. It is interesting 
to note that this spectral index is slightly softer than the 2.45 ± 0.07 we obtained for the 
low flux (17 hours of observation, ~0.4 c.u.) in the energy range 0.1-6 TeV. This spectral 
shape difference might be caused by a possible softening of the spectra above 1-2 TeV. Note 
that the fits to the spectra measured by MAGIC are mostly constrained by the points below 
2 TeV (due to low photon statistics at the highest energies), while the fits to the spectra 
measured by HEGRA are mostly constrained by the points above 1 TeV. Certainly, a factor 
contributing to the softening of the spectra at the highest energies is the 7-ray extinction due 
to pair productio n by in teraction with the Extragalactic Background Light (EBL). According 



to 



Kneiske et al I (120041 ). the attenuation of 7-rays coming from Mrk 501 is ~30% at 1 TeV 



and ~50% at 10 TeV, while it is <15% for energies below 0.5 TeV. 



An exponential cutoff of ^^5 TeV was suggested in lAharonian et al I (j200l[ ). although the experimental 
data are perfectly compatible with both hypotheses, the simple PL {y^ /NDF = 12.9/14; P = 0.53), and the 
PL with exponential cutoff {x^/NDF = 9.2/13; P = 76). 
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4.3. Spectra during active nights 

The differential photon spectra from June 30 and July 9 were fitted with the simple PL 
described in eq. |5] as well as with the log-parabolic function: 

dF f E \ -a-b-fo3io( „ 3g^y ) 

dE = ^'\-^2T^) ■ 

Here, Kq is a normalization factor, a is the spectral index at 0.3 TeV, and 6 is a curvature 
parameter (for 6>0 the spectrum hardens/softens at energies below/above 0.3 TeV) . The log- 



parabola is a simple function to describe curved spectra and, as pointed out by lMassaro et al 



(12004 . 120061 ) . can be directly related to the intrinsic physical processes occurring in the source. 
The results of both, the power law and the log-parabolic functions are reported in Table 
Since the log-parabolic fit describes the data more accurately than the PL fit, this suggests 
that the differential photon spectra of the two flaring nights are curved. The spectra and 
their flts are shown in Fig. [TH Note that the peak of the energy spectrum is located in the 
covered energy range in both nights. 

The peak location in a spectrum described by eq. E] is given by: 

Epeak = 10^ -0.3 TeV (7) 

with an associated uncertainty of 



InlO / /a — 2\2 ~ 2\ 

^Epeak = Epeak " ' V Vaa + Hfe " ^ J " 2\4fe " ^ J (8) 

where Vaa, Vhb and Vab are the coefficients of the covariance matrix. Using these equations 
one finds that the peak locations are 0.43±0.06 TeV and 0.25±0.07 TeV for the spectra 
measured during June 30 and July 9, respectively. It should be noted that these spectra 
are not corrected for EBL absorption, and are therefore not intrinsic to Mrk 501. After 
correction for EBL absorption the spectral peaks are shifted towards higher energies (see 

Sect. ra . 

The detection of a spectral curvature and m easurement of the peak location was first 



reported by the Whippl e and CAT collaborations (ISamuelson et al Ill998l : iDjanati-Ataiet al 



19991 : iPiron Il200d . l200ll ) based on 1997 Mrk 501 data. However, in those studies the spectra 
were not corrected for the EBL absorption. Such correction is relevant because the measured 
curvature in those energy spectra occurred essentially above 1-2 TeV, where current EBL 
models predict a >40% attenuation. Below 1 TeV, little (if at all) curvature could be seen 
within the quoted la statistical errors. Therefore, we suggest that the curvature and peak 
location reported in those studies may be significantly affected (if not dominated) by the 
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EBL absorption. On the other hand, the spectral curvature in the MAGIC data is dominated 
by points below 1 TeV, since the higher-energy points hardly constrain the fit due to their 
substantially larger error bars. The curvature we measure thus is affected, but not dominated 
by EBL absorption. 

4.4. Intra-night spectral variations 

As discussed in section 13. 3[ during the active nights of June 30 and July 9 the VHE 
emission of Mrk 501 can be divided into a 'stable' (pre-burst) and 'variable' (in-burst) part. 
In order to study potential changes in the spectral shape, we derived the differential photon 
spectrum for the two parts of each night (Fig. [5]). Since each spectrum is based on <l/2 hr 
exposure this procedure certainly increases the statistical errors. The four spectra were fitted 
with the log-parabolic function in eq. [6l which is preferred over the simple PL function (see 
sect. 14.31) . The results of the fit, as well as other relevant information (e.g., net observing 
time, significance of the signal, goodness of fit) are reported in Tabled The spectra and the 
corresponding fits are plotted in Fig. [151 In both nights, there is marginal ( la) evidence 
for a spectral hardening during the flare. 

We also studied the time-evolution of the hardness ratio, defined as the ratio 
F(1.2-10 TeV)/F(0.25-1.2 TeV) and which is computed directly from the LCs shown in Fig. 
|6] and Fig. [71 The resulting graph is shown in Fig. [161 The hardness ratios for the pre-burst 
and in-burst part are quantified by means of a constant fit. In both nights the hardness ratio 
is somewhat larger (l-2cr) in the in-burst than in the pre-burst part, in agreement with the 
observed spectral hardening (see Fig. [TB]) . It is worth noting that the hardness ratios for 
the pre-burst and in-burst time windows of June 30 are statistically compatible with being 
constant, those of July 9 are much less so, as shown in the insets of Fig. [161 

The evolution of the hardness ratio with the emitted flux above 0.25 TeV is shown in 
Fig. [IZl Both nights show some evidence for a larger spread in the in-burst part than in 
the pre-burst part. The evolution of the in-burst points from June 30, however, is somewhat 
chaotic, while the evolution of the in-burst points from July 9 shows a clear loop pattern 
rotating counterclockwise. The physical interpretation of this feature is given in section 15.31 
Concluding, a spectral hardening with increased emission characterizes the VHE emission of 
Mrk 501 also at short timescales. 
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5. Discussion 

In this section we discuss the VHE light curve of Mrk 501 for May - July 2005 in 
comparison to previous lACT observations made from the years 1997 through 2000. We 
also discuss the observed rapid flux variabihty in the framework of a basic model of gradual 
electron acceleration. The broadband spectral features of Mrk 501 and the intrinsic VHE 
spectra from different activity states are discussed in the framework of a basic SSC model. 
We also model the ensuing intrinsic VHE spectra - corresponding to different activity states 
- in an SSC framework. 



5.1. Historical light curve 



It is interesting to examine Mrk 501's VHE activity in 2005, as measured by MAGIC in 



persp ective to that recorded in 1997-2000 by other lACTs s uch as HE GRA CTl (jKranich 



2000), HE GRA CT Svstem flMulti-Messenger-Group I I2006D. Whipple (IQuinn et al 



19991 



and CAT JPianati-Ataiet al Ill999l :" iPiron II2OOII ). The long-term Mrk 501 lightcurve covering 



the years 1997 to 2005 is shown in Fig. O For easier comparison among instruments 
covering different energy ranges, the integral flux values are given in Crab Units (c.u.). 
The mean fluxes for each year and instrument are shown in the insets of Fig. [181 They 
were derived using the exposure times as statistical weights (whenever this information was 
available), and excluding data points > 3a away from each mean value, to permit a better 
comparisorj^. For MAGIC, mean fluxes were computed in the other lACTs' energy bands. 
In 1997, when Mrk 501 was much brighter than the Crab, Whipple and CAT fluxes were in 
mutual agreement, but they significantly differed from the HEGRA CT System and HEGRA 
CTl fluxes; the reasons being (i) the highest-flux nights (May- July 1997) happened to be 
covered by HEGRA CT and CTl but not by Whipple and CAT; and (ii) the highest threshold 
energies of HEGRA CT and CTl, together with the fact that when Mrk 501 is active its 
spectrum becomes significantly harder than the Crab's (see Sects. 14. H and l4.2p . hence fluxes 
measured in Crab units are larger at higher energies. In 1998 and 1999 the various lACT 
data agreed on a flux of ~0. 15-0. 20 c.u., i.e. an order of magnitude lower than in 1997. In 
2000, March through May, the mean VHE flux of Mrk 501 increased to a level of 0.35±0.09 
(HEGRA CTl) and 1.19±0.17 (Whipple), the main reason for the discrepancy being that 
HEGRA CTl missed the nights with highest flux measured by Whipple. In 2005, May 



Some HEGRA CTl, Whipple and CAT fluxes are negative, unsurprisingly because many of the fluxes 
from these instruments are just ^ 2a measurements. For the HEGRA CT System we have no information 
about negative values, hence the corresponding mean fluxes may be slightly overestimated. 
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through July, the mean basehne VHE flux from Mrk 501 was ~0.5 c.u., significantly lower 
than in 1997 and 2000, but higher than in 1998 and 1999. 



A 23-day fl ux periodicity was claimed by iKranich I (120001 ) using data from HEGRA 
CTl. Recently, lOsone I (120061 ) confirmed the 23-day periodicity at VHE frequencies and 
extended it to X-rays (based on ASM/RXTE data) using a more sophisticated time analysis. 
This periodic modulation in the VHE emission of Mrk 501 may be evidence of a binary 
black hole system with separation of the order of the gravitational radius, as suggested by 
Rieger fc Mannheim I (l2000l . l2002l ). or selective absorption of 7-rays in the radiation of a hot 



spot or biting in the inner part of the accretion disk, as pointed out by lBednarek fc Protheore 



( ll997bl ). This periodicity was not seen in the 1998-2000 campaigns, when the source was 
apparently not very active. This might indicate that such a periodicity in the emitted flux 
occurs only when the source shows very high activity. The MAGIC 2005 data did not have 
the required coverage for such a timing analysis. 



5.2. Interpretation of the measured fast flux variations and energy dependent 

time delays 

The very short flux-doubling time and the energy-dependent time delays of the flux 
variations can give us information on the acceleration processes occurring in Mk501. In this 
section we argue that gradual electron acceleration in the emitting plasma can provide a 
natural explanation of the observed time structures. 

i) Flare decay timescale. 

Let us assume that the maximum energies of electrons accelerated in the relativistic 
blob are determined by their radiation energy losses on synchrotron and IC processes, as 
expected in the SSC model. The acceleration time should then equal the energy loss time 
scale 

''"acc ^cool; (9) 

where 

T^,, = E'JP^,, = E'J{icE'jR^)^^.lE'J^B s. (10) 

Here E'^ is the electron energy (in TeV) in the blob. Pace is the rate of energy gain during 
the acceleration process, B (in G) is the magnetic field strength in the acceleration region, 
i?L is the electron Larmor radius, and ^ is the acceleration efficiency. The cooling time of 
electrons can be expressed by 
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111 
= — + — . (11) 

''"cool ^syn ''"IC 

Then, the cooling time can be expressed by only the synchrotron cooling time 

Tcool = Tsyn/(1 +^?), (12) 

where 77 = Tgyn/ric, Tgyn = E'^/Pgy^, P^yn = 47rcrTcf/B7e/3, ctt is the Thomson cross section, 
c the velocity of light, Ub = B'^/Stt is the energy density of the magnetic field, 7e = E'^/mc^, 
and m is the electron rest mass. The parameter t] corresponds to the ratio of the power 
emitted by electrons in IC and synchrotron processes, respectively (which is reciprocal to 
the coresponding cooling times Tsyn/ric). The MAGIC 7-ray and corresponding RXTE/ASM 
X-ray data permit to constrain to < 0.7. The modelling of the SED presented in section [531 
(see. Fig. [21]), however, suggests that t] is more likely of the order of ~ 0.2. We therefore used 
7^=0.2 in all the following estimates. By comparing Eqs. [9|T0l and fT2l and setting E'^ ^ Etby/S 
TeV (where 6 is the Doppler factor of the relativistically-moving emitting plasma blob and 
ii^TeV is the electron energy (in TeV) in the observer's frame, we obtain the condition on the 
acceleration efficiency of electrons, 

^ ^ 10-'BEly{l + r])S-'. (13) 

If the observed decay timescale of the flare, Tf, is also due to radiative processes, then, 

Tcool = Ti5 . (14) 

Inserting Eq. [T2] into Eq. [TH and expanding Tgyn and -Pgyn, permits to estimate B in the 
cooling region: 

B ^ 11.2[(l + r7)-(ETevT-f)]-i/2 G, (15) 
with Tf in seconds. If electron acceleration and cooling are co-spatial, Eq. [13] becomes 

^ ^ 1.1 X 10-2 ^TeV U + ^J (Ig) 

For the parameters of the July 9 flare MrkSOl, the characteristic flux variability time Tf 3 
minutes (see Sect. 13.31) and i?TeV ~ 1, we obtain E ^ 0.8 G and ^ fa 0.9 x lO"^^"^. Note, 
that this estimate of the magnetic field in the emission (acceleration) region of Mrk 501 is 
consistent with previous estimates, based on a homogeneous SSC model and the assumption 
of a very short variability time scale such as for the April 15-16 1997 fiare (see e.g.. Fig 3c in 
Bednarek fc Protheore I Jl999h l 
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ii) Energy- dependent time delay in peak flare emission. 

The time delay between the peaks of F(<0.25 TeV) and F(>1.2 TeV) during the July 
9 flare, can be interpreted as due to the gradual acceleration of electrons in the relativistic 
blob. As reported in section 13. 3[ under the assumption that the shape of the flares is the 
same in the two energy ranges, the time delay is Ath-l = 4 ± 1 minutes. 

Within the above framework, the time delay should correspond to the difference between 
the acceleration times of electrons to energies -ExeV ~ 0-25 and -ExeV ~ l-^) according to 

^a"cc-^acc = A^H-L^- (17) 

Assuming for simplicity that the electron energies are determined by the energies of the 
emitted VHE radiation, i.e. E'^RiE'^^y/S and E'J^^E'^y/5, we can use Eq. [10] to model the 
acceleration time of electrons in eqJTTl Finally, by reversing Eq. [T71 we get another limit on 
the electron acceleration efficiency, 



L 



Applying the observed values of -ExeV; -^TeV; ^''"h-l and the estimate of B from Eq. [TS] (which 
is valid if electron acceleration and cooling are co-spatial) we obtain x 10~'^5~^, which 

is roughly the same value as determined above using the exponential flux decay time. Note 
that diffusion and/or spatial dishomogeneities might increase the volume where electrons 
cool down, thus making the magnetic fleld in the cooling region lower than that of the 
acceleration region. That would imply that the estimated values for the parameter ^ using 
equations [16] and [18] are upper and lower limits, respectively. 

We conclude that the time delay of the flare peak emission in different ranges of energy 
can result from the gradual acceleration of t he emitting electro n s in t he blob. The inferred 
blazar Doppler factors, 5 ~ 10 — 15 (e.g. Costamante et al ( 20021 )). imply a relatively 



inefficient acceleration, ^ ~ 10^^. This value is significantly lower than required by the 
observations of 7-ray emission from the pulsar wind nebulae in which leptons are accelerated 
in the shock wave in the relativistic pulsar wind. For example, in the Crab Nebula ^ has to be 
of the order of ~ 0.1, since leptons are accelerated clearly above 10^ TeV (approximately 10% 
of the maximum available potential drop through the pulsar magnetosphere). Therefore, the 
physics of the acceleration process in the relativistic jets of BL Lacs and relativistic shocks 
in the pulsar wind nebulae may differ significantly. 

A somewhat more speculative issue that blazar emission permits to explore concerns 
non-conventional physics. Energy-dependent arrival times are predicted by several models 
of Quantum Gravity, which quantify the first-order effects of the violation of Lorentz sym- 
metry. One could, therefore, speculate that the observed time difference is explained by such 



-24- 



models, although source-inherent effects could certainly not be excluded. A more detailed 
investigation of such interpretations of our data is still going on. 



5.3. Interpretation of the spectral shape variations 



The observed correlation between s pectral shape and ( bolometric) luminosity is natu- 
rally accounted for in the SSC scenario. iPian et al I (119981 ): iTavecchio et al I (120011 ) discuss 
the SED variations of Mrk 501 during the giant 1997 flare. During the 1997 flare the 0.1- 
200 keV band synchrotron spectrum became exceptionally flat (photon spectral index a<l), 
peaking at <^ 100 keV - a shift to higher frequencies by a factor of 100 from previous, more 
quiescent states. The VHE data (from the Whipple, HEGRA, and CAT telescopes) showed 
a progressive hardening from the baseline state (a>2) through a more active state to a flar- 
ing state (a~2). In the SSC scenario, these flux- dependent spectral changes implied that 
a drastic change in the electron spectrum caused the increase in emitted power: a freshly 
injected electron population has a flatter high-energy slope and a higher maximum energy 
than an aging population, which cause a shift of the SED to higher frequencies. 



In section U3] we reported that the spectrum of Mrk 501 not only hardens on long time 
scales, with the overall emitted flux, but also during the shorter events of June 30 and July 9. 
The burst from July 9 showed a remarkable variability, and the evolution of the hardness ratio 
with the flux (right-hand plot of Fig [T71) contains valuable information about the dynamics 
of the source. In the pre-burst phase, the hardness ratio does not vary significantly; yet 
during the burs t phase, it varies following a clear loop pattern rotating counterclockwise. As 
pointed out by iKirk fc Mastichiadis I (119991 ). one expects to have this behaviour for a flare 
where the variability, acceleration and cooling timescales are similar; which implies that, 
during this flare, the dynamics of the system is dominated by the acceleration processes, 
rather than by the cooling processes. Consequently, the emission propagates from lower to 
higher energy, so the lower energy photons lead the higher energy ones (that is the so-called 
hard lag). This indeed agrees well with the argumentation given in section [3721 where the 
time delay between -E>i.2 TeV and -E<o.25 Tev is shown to be consistent with the gradual 
acceleration of the electrons. 



In a systematic study performed by iGliozzi et al I (120061 ) using X-ray data from 1998 
to 2004, this behaviour was not observed on more typical (longer) flares, where actually the 
opposite behaviour (clockwise rotation) was indicated. This might point to the fact that 
these physical processes might be responsible only for the shortest flux variations, and not 
for the variability on longer timescales. 
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5.4. Interpretation of the increased variability with energy 

In the SSC framework, the variabihty observed in the VHE emission brings informa- 
tion about the dynamics of the underlying population of relativistic electrons (and possibly 
positrons). In this context, the general variability trend reported in section [3^ is interpreted 
by the fact that the VHE 7-rays (as well as the X-rays) are produced by more energetic par- 
ticles, that are characterized by shorter cooling timescales; causing the higher variability 
amplitude observed at the highest energies. It is worth noticing that such an injection of 
high energy particles would produce a shift in the IC peak, which is indeed observed during 
this observing campaign, as reported in sections HI and 1^31 



5.5. Spectra corrected for the EBL absorption 

In this section we correct the measured spectra for the EBL absorption. To this purpose 
we use Kneiske et al. (2004: 'Low'). A correction using the EBL models Aharonian et al. 
(2006) and and Primack et al. (2005) gave very similar results, while using Kneiske et al. 
(2004: 'Best') and Stecker et al. (2006,2007) gave slightly larger energy fluxes above 1 TeV0. 
Fig. [T^ shows the spectra from the active nights June 30 and July 9 before/after correction 
for the EBL absorption. In spite of our proximity to Mrk501, the effect of the EBL is not 
negligible, and the spectral peak moves to higher energies. 

The location of the spectral peaks (calculated using eqs. [7]and[HD are shown as a function 
of F (>0.15 TeV) in Fig. [20] for the flaring nights before and after the EBL correction. The 
flgure seems to indicate a displacement of the peak location with the increasing flux, yet the 
error bars are too large to be conclusive. On the other hand the peak location is certainly 
at <0.1 TeV when Mrk 501 is in a low state (see sect. 14.21) . Hence, there is evidence for an 
overall peak location versus luminosity trend. 

Fig. [21] shows the June 30 and July 9 spectra, as well as the mean 'high', 'medium', and 



low '-fl ux spectra (see sect. 14. 2p . EBL-de-absorbed using the 'Low' model of [Kneiske et al 



( I2OO4J ). The X-ray fluxes measured by RXTE/ ASM, and the optical flux observed by 



the KVA Te l escope , are also shown. The optical flux from the host galaxy, estimated by 



Nilsson et all (|2007[ ) to be 12.0 ±1.2 mJy, has been subtracted. 



The best-fitting one-zone, homogeneous SSC models of Mrk50rs intrinsic spectrum for 
the highest state of the source (corresponding to the active night of June 30) and for the 



The difference increases witli tfie energy, being ~50% at 5 TeV. 
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lowest state (the mean spectrum corresponding to the 'low'-flux bin, see Fig. [T3|) are dis- 
played in Fig. [211 The fit parameters (electron population's break and max/min energies, 
high/low-i? spectral slopes, normalization: _Ebr, -Emax, -Emin, ^i, n2, K; plasma blob's ra- 
dius, magnetic field, Doppler factor R, B, 6) are reported in Table M (see Tavecchio et al. 
2001 for details on the model). It should be noted that two different fits to the high-state 
spectrum are possible - whose main differing parameters are, respectively, S,B=25, 0.23 G 
(solid black curve) and 50, 0.053 G (dashed black curve) - which show that fairly different 
synchrotron peaks are possible within our X-ray and (EBL-corrected) TeV data. Spectrally 
more extended X-ray data would probably have solved the degeneracy. The optical data, 
too, do not lift the degeneracy: once the optical light contribution of the underlying host 
galaxy is subtracted, the observed energy flux is rather compatible with the SED models for 
the different activity state^. The flt to the 'low'-state spectrum is characterized, perhaps 
unsurprisingly, by a change of the internal physical conditions of the emitting plasma blob 
rather than by a change of its bulk attributes (blob size and relativistic Doppler factor): the 
low state is characterized, with respect to the high state, by lower max/break energies and 
normalization of the electron population and by a somewhat stronger magnetic fleld. One 
nice consistency feature of all the flts is that, in all cases, the radius of the plasma blob, 
R=10^^ cm, implies a crossing time tcr=R/Sc, comparable to that inferred from the observed 
duration (~20 minutes) of the flare, Atflare- 

The SED models for Mrk 501 derived and discussed in this section can be compared 
with some previous published models, like Plan et al. (1998): 5, i?~15,0.8G from one-zone 
SSC modeling of 1997 SEDs, with flatter /steeper electron distributions for active/ quiescent 
phases; Tavecchio et al. (1998): 6,Bc^8 — 20,0.5 — IG from one-zone SSC modeling of 
historical quiescent SED and 6, -B~7, IG (and very high break energy) for the active SED; 
Bednarek & Protheroe (1999): 5, S~12 - 36,0.07 - 0.6G from 1997 SEDs modeled with 
one-zone SSC requiring 7-ray transparency of the emitting blob; Kataoka et al. (1999): 
S, i?~15, 0.2G from one- zone SSC modeling of simultaneous 1996 SED; Katarzyhski et al. 
(2001): 6,B = 14, 0.2G from SSC modeling of non-simultaneous broadband SED. The SED 
of Mrk 501 can also be modelled by less conv entional approaches, requiring magnetic flelds 



in the emission region smaller than 0.005G, iKrawczynski I (120071 ). We should, however 



remember that for this highly variable source, constraints derived for some epochs may not 
apply to other epochs (the simple one-zone model has 9 free parameters!): in particular, 
most published models refer to the giant 1997 flare hence comparisons with our results may 
not be straightforward. However, Tavecchio et al. (2001) modelled different emission states 



The measured optical flux might have contributions from regions outside the one producing the radiation 
at X-ray and 7-ray frequencies. 
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of MrkSOl in 1997, 1998, and 1999 by just changing the electron energy distribution (slopes, 
break energy, number density) and keeping the other parameters frozen, similarly to what 
was done here (see in Tabled the t wo states represented by so lid lines in Fig. [51]). It is also 
worth mentioning the work done by iKrawczynski et al I (120021 ) . in which the energy spectra 
(for different days) were modelled using a time- dependent code. 

It should be remarked that in TeV blazars, while the bright and rapidly variable VHE 
emission implies that at the scales where this emission originates (~0.1 pc from the jet apex) 
the jet is highly relativistic (5~10-20 with no EBL-absorption correction, and 50 with 
correction), at VLBI (~1 pc) scales the jets are relatively slow (see Ghisellini et al. 2005, 
and references therein). Hence, to reconcile the high 5's derived from VHE data with the 
much lower (5's derived from VLBI radio measurements, the jets of TeV blazars must either 
undergo severe deceleration (Georganopoulos & Kazanas 2003) or be structured radially as 
a two- velocity, inner spine plus outer layer, flow (Ghisellini et al. 2005). 



6. Concluding remarks 

In this work, we have undertaken a systematic study of the temporal and spectral 
variability of the nearby blazar Mrk 501 with the MAGIC telescope at energies > 0.1 TeV. 
During 24 observing nights between May and July 2005, all of which yielded significant 
detections, we measured fluxes and spectra at levels of baseline activity ranging from <0.5 
to >1 C.U.. During two nights, on June 30 and July 9, Mrk 501 underwent into a clearly 
active state with a 7-ray emission >3 c.u., and fiux-doubling times of ~2 minutes. The ~20- 
minute long flare of July 9 showed an indication of a 4±1 min time delay between the peaks 
of F(<0.25 TeV) and F(>1.2 TeV), which may indicate a progressive acceleration of electrons 
in the emitting plasma blob. An overall trend of harder spectra for higher flux is clearly seen 
on intra-night, night-by- night, and longer-term timescales. The VHE 7-ray variability was 
found to increase with energy, regardless whether the source is in active or quiescent state, 
and it is significantly higher than the variability at X-ray frequencies. A spectral peak, at a 
location dependent on source luminosity, was clearly observed during the active states. All 
these features are naturally expected in synchro-self- Compton (SSC) models of blazar VHE 
emission. There are no simultaneous good quality X-ray measurements during the MAGIC 
observations. As a consequence, the SSC model of the X-ray/VHE SED of Mrk 501 in an 
active state is not unequivocally constrained, but it still restricts the emitting plasma blob 
to have Doppler factors in range of 25-50 and magnetic fields in the range 0.05-0.5 G. 
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Table 1. MAGIC observation of Mrk 501 



MJD 
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Tabs ^ 
(11,) 




COTTIO 




ph\ 
s ' 


^>0.15 TeV 

(Crab Units) 
















P h 


53518.980 


0.75 


19.10-28.95 


6.44 


1.19 


± 


0.25 


0.37 


± 


0.08 


2.63 ± 


0.48 


2, 


.17 


± 


0. 


,25 


2.7/8 


95.2 


53521.966 


1.85 


9.97-30.10 


8.90 


1.51 


± 


0.17 


0.47 


± 


0.05 


2.94 ± 


0.33 


2, 


.61 


± 


0, 


16 


10.8/7 


15.0 


53524.969 


0.58 


19.18-27.73 


6.98 


2.04 


± 


0.29 


0.64 


± 


0.09 


3.71 ± 


0.53 


2, 


.47 


± 


0, 


,23 


1.6/6 


95.0 


53526.975 


0.98 


9.96-28.94 


8.69 


1.63 


± 


0.22 


0.51 


± 


0.07 


3.26 ± 


0.38 


2, 


.49 


± 


0, 


17 


3.8/9 


92.4 


53530.973 


0.47 


15.22-22.32 


6.52 


1.53 


± 


0.32 


0.48 


± 


0.10 


2.28 ± 


0.65 


1, 


.97 


± 


0. 


,49 


1.1/3 


78.9 


53531.959 
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15.21-25.15 


6.98 


1.29 


± 
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± 
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2.69 ± 
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0. 


,30 
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2.41 ± 
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.34 


± 


0. 
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± 
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2.46 ± 
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0, 


,19 
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0.45 
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± 
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44.8 
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3.43 


9.95-30.07 
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2.69 
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0.85 
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4.45 ± 


0.24 


2, 


.42 


± 
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,06 
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45.3 
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2.68 


9.95-29.93 


7.01 
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0.24 
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1.36 ± 
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± 


0, 
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1.25 


± 


0.10 


0.39 


± 


0.03 


2.08 ± 


0.19 


2, 


.46 


± 


0. 


,14 


8.2/8 


41.4 


53548.931 


0.87 


9.98-20.68 


6.12 


1.21 


± 


0.25 


0.38 


± 


0.08 


2.39 ± 


0.38 


2, 


.28 


± 


0. 


,27 


0.6/6 


99.6 


53551.905 


1.09 


12.86-25.15 


32.02 


11.08 ± 0.32 


3.48 


± 


0.10 


17.37 ± 0.51 


2, 


.09 


± 


0. 


,03 


26.2/11 


0.6 


53554.906 


0.68 


15.21-22.32 


12.52 


3.52 


± 


0.30 


1.11 


± 


0.09 


5.91 ± 0.47 


2, 


.26 


± 


0. 


,11 


3.9/9 


92.1 


53555.914 


0.44 


12.85-22.32 


6.08 


1.27 


± 


0.34 


0.40 




0.11 


2.96 ± 0.62 


1 


.97 




0, 


,29 


1.9/6 


92.5 


53557.916 


0.54 


12.84-19.06 


8.40 


2.25 


± 


0.32 


0.71 


± 


0.10 


3.91 ± 0.48 


2, 


.30 


± 


0, 


,21 


6.5/7 


48.5 


53559.920 


0.98 


9.94^17.22 


10.05 


1.85 


± 


0.23 


0.58 


± 


0.07 


3.10 ± 0.33 


2, 


.25 


± 


0, 


,13 


8.4/8 


39.9 


53560.906 


0.76 


9.96-19.07 


24.39 


9.93 


± 


0.38 


3.12 


± 


0.12 


14.35 ± 0.56 


2, 


.20 


± 


0. 


,04 


22.5/11 


2.1 


53562.911 


1.63 


9.94-16.79 


11.08 


2.19 


± 


0.37 


0.69 


± 


0.12 


2.83 ± 0.30 


2, 


.34 


± 


0. 


,13 


14.1/8 


8.2 


53563.921 


0.85 


9.94-15.16 


18.69 


5.53 


± 


0.28 


1.74 


± 


0.09 


7.89 ± 0.39 


2, 


.25 


± 


0. 


,06 


11.5/9 


24.3 


53564.917 


0.34 


9.94-15.18 


8.91 


2.89 


± 


0.46 


0.91 


± 


0.15 


4.88 ± 0.56 


2, 


.27 


± 


0, 


,20 


5.4/6 


49.7 


53565.920 


2.57 


9.95-28.93 


11.62 


1.71 


± 


0.13 


0.54 


± 


0.04 


2.73 ± 0.22 


2, 


.49 


± 


0, 


12 


10.7/8 


21.6 


53566.953 


1.91 


9.99-30.10 


11.63 


1.33 


± 


0.11 


0.42 


± 


0.04 


2.16 ± 0.20 


2, 


.28 


± 


0. 


,13 


7.4/10 


69.0 



*Net observation time after removing bad-quality runs. 

''Zenith Angle range covered during the observation. 

'^Combined significance of detected signal in the 0.1-10 TeV band. 

■^Integrated flux above 0.15 TcV. 

"Normalization factor of the PL fit. 

''Slope of PL fit. 

value and number of degrees of freedom of the power-law fit. 
''Chance probability for larger values. 
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Table 2. Flare model parameters: integral emission above 0.15 TeV. 



Date 


oos 


'-'comb 


a'' 


a 


b 


c 


d 




p e 




(h) 


(sigma) 


^ cm^ -s ' 


(Crab Units) 


^ crri^ • s ^ 


is) 


is) 




{%) 


June 30 


0.63 


24.7 


10.80±0.48 


3.39±0.15 


13.2±4.7 


81±41 


50±23 


20.0/15 


17.3' 


July 9 


0.36 


19.6 


7.39±0.48 


2.32±0.15 


20.3±3.3 


95±24 


185dz40 


4.2/7 


75.8 



''Net observation time during variable emission (right part of the graphs). 

''Combined signal significance from variable emission (right part of the graphs) in 0.1-10 TeV band. 
'^Integrated flux above 0.15 TeV for the steady emission (left part of the graphs). 
'^X^ value and number of degrees of freedom of the fit with eq. |2] 
"Chance probability of having larger values. 

'if the points after 22:44 are not taken into account, the coefficients from the fit remain the same, and the probability increases 
up to 52.7% ix^/NDF = 11.0/12). 



Table 3. Flare model parameters for June 30: differential emission. 



Energy Range a'^ 

(TeV) (1^^) 


a b 
(Crab Units) ('"'" P'^) 


c 

is) 


d 

is) 


X^/NDF'^ 


P <= 
{%) 


0.25-0.6 3.30±0.23 


3.0±0.2 7.5±2.8 


110±57 


61±26 


5.2/6 


51.8 



'^Integrated steady emission flux (left part of the graphs) in specified energy range. 
'^X^ value and number of degrees of freedom of the fit with eq. |2] 
''Chance probability of having larger x^ values. 



Table 4. Flare model parameters for July 9 resulting from a combined fit to all LCs from 
Fig. El using equation [2] with c=d. The overall x^/NDF = 14.0/12 (P=0.3) 



Energy Range 
(TeV) 


a'' 

('°"'° "") 

^ cm'^-s ' 


a 

(Crab Units) 


b 


c 

{^) 


, ,hC E 0.15-0.25TeVb 
1^0 ~ t^o 

is) 


0.15-0.25 


4.23±0.49 


2.48±0.28 


8.6±3.7 


143±92 


± 68 


0.25-0.6 


2.55±0.24 


2.32±0.09 


9.3±2.5 


95±28 


7 ± 36 


0.6-1.2 


0.53±0.10 


1.96±0.37 


2.7±0.9 


146±56 


111 ± 91 


1.2-10 


0.23±0.06 


1.51 ±0.39 


4.0±0.9 


103±19 


239 ± 40 



''Integrated steady emission fiux (left part of the graphs) in specified energy range. 

^t^^ ^ 0.i5-0.2jTeV ^j^g £^j. £j^g ^j^g energy range 0.15-0.25 TeV. This is used as a reference value, and the error 

of this quantity is not taken into account. 
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Table 5. Flare model parameters for July 9 resulting from a combined fit to all LCs from 
Fig. [8] using equation [2] with c = d, and with a common to for all LCs. The overall 

X^/NDF = 25.6/15 (P=0.04) 



Energy Range 
(TeV) 




a 

(Crab Units) 


b 

,10^^° ph-. 
^ cm-^-s ' 


c 

(s) 


, ,LC E 0.15-0.25TeVb 
(s) 


0.15-0.25 


4.23±0.49 


2.48±0.28 


5.4±2.2 


301±210 


± 42 


0.25-0.6 


2.55±0.24 


2.32±0.09 


5.7±1.5 


162±63 


± 42 


0.6-1.2 


0.53±0.10 


1.96±0.37 


2.6±0.8 


153±56 


± 42 


1.2-10 


0.23±0.06 


1.51 ±0.39 


3.9±1.0 


97±22 


± 42 



'^Integrated steady emission flux in specified energy range (left part of the graphs). 

^^0^ ^ 0.i5-0.25Tey ^j^^ ^^^^ energy range 0.15-0.25 TeV. This is used as a reference value, and the error 

of this quantity is not taken into account. 



Table 6. Stacked analysis: mean spectral parameters. 



Flux Level"" 


^ ohs 
(h) 


(deg) 


o c 
-J comb 

sigma 


-r>0.15 TeV 
^ cm-^-s ^ 


F>0.15 TeV 
(Crab Units) 


, 10~^° ph ^ 
em'' -s-O.iTeV ' 






ph 

(%) 


Low 


17.2 


9.96-30.1 


16.7 


1.24± 0.08 


0.39 ± 0.02 


2.31±0.13 


2.45 ± 0.07 


7.8/7 


34.6 


Medium 


11.0 


9.95-30.0 


22.8 


2.11± 0.09 


0.66 ± 0.03 


3.57±0.15 


2.43 ± 0.05 


2.9/7 


89.4 


High 


1.52 


9.95-22.3 


21.7 


4.62± 0.21 


1.45 ± 0.07 


7.13±0.32 


2.28 ± 0.05 


4.8/7 


68.7 



™See section [4.21 for definition of fiux levels. 

''Net observation time after removing bad quality runs. 

''Zenith Angle range covered during the observation. 

'^Combined significance of detected signal in energy range 0.10-TeV-lO TeV. 

Integrated fiux above 0.15 TeV. 
''Normalization factor of the power-law fit. 
^Slope of the power-law fit. 

value and number of degrees of freedom of the power-law fit. 
'"Chance probability of having values. 



Table 7. PL (eq. [5]) and Log-parabolic (eq. [6]) fit results for June 30 and July 9: 

night-integrated spectra. 



Fit performed with eq. \5\ Fit performed with eq. \E\ 

Date K^) a x^/NDF P ^ Wo a h x^/NDF P 

, 10-1" ph . f 10-"' ph ,0/-) 
^cm'^-s-O.STeV ' ^ cm'-' ■ s-0.3TeV ' Y.' 

June 30 17.4±0.05 2.09±0.03 26.1/11 0.6 18.6±0.06 1.89±0.06 0.35±0.09 6.1/10 80.1 

July 9 14.3ili0.06 2.20ili0.04 22.5/11 2.1 15.5ifc0.07 2.06ifc0.07 0.36ifc0.16 15.2/10 12.5 



'^Chance probabiUty for larger values. 



Table 8. Log-parabolic (eq. [6]) fit results for June 30 and July 9, 2005: pre-burst and 

burst spectra. 



MJD Tois " Sccnb'' i^>0.15 Tev"" i^>0.15 TeV Kq a b x'iNDF P ^ 

Start ih) Sigma {^^^^) (Crab Umts) { JCLt.v ) 

53551.905 0.46 22.3 10.99±0.48 3.46±0.15 19.8±1.0 1.97±0.08 0.27±0.14 8.2/9 51.2 
53551.924 0.63 24.7 11.15±0.43 3.50ili0.14 17.2ili0.8 1.87ili0.08 0.34±0.13 13.8/10 18 .1 

53560.906 040 15^2 7.64 ± 0.48 2.40 ± 0.15 12.7±1.1 2.11±0.12 0.57±0.34 6.4/8 59.8 
53560.923 0.36 19.6 12.39 ili 0.60 3.89 ± 0.19 19.3ilil.3 2.00ili0.10 0.44ifc0.23 8.9/8 35 .2 

''Net observation time after removing bad-quality runs. 
'^Combined significance of detected signal in energy range 0.1-10 TeV. 
'^Integrated 7-ray flux above 0.15 TeV. 
'^Chance probability for larger values. 



Table 9. SSC model parameters 



spectrum 


O'min 


7br 




nl 


n2 


B 


K 


R 


Doppler factor 














Gauss 


particle/cm^ 


cm 




June 30 


1 


10'^ 


10^ 


2 


3.5 


0.23 


7.5- 10"* 


10" 


25 


June 30 (bis) 


1 


5 ■ 10^ 


10^ 


2 


3.5 


0.053 


7.0- 10'' 


1015 


50 


Low flux 


1 


10^ 


5 ■ 10^ 


2 


3.2 


0.31 


4.3- 10" 


1015 


25 
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Fig. 1.— Multi-frequency LC during the MAGIC observations of Mrk 501 (May- July 2005). 
Top) MAGIC flux above 0.15 TeV. The Crab flux is also shown for comparison (lilac dotted 
horizontal line). Middle) RXTE/ ASM 2-10 keV flux. Bottom) KVA ~1.5-2.5 eV flux. Error 
bars denote la statistical uncertainties. The X-ray/optical data were selected to match the 
MAGIC data within a time window of 0.2 days. 
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Fig. 2. — VHE versus X-ray (left) and optical (right) flux correlation during the MAGIC 
observation campaign. The data points are the same as in Fig. [H The blue dotted lines 
denote constant fits; black solid and green dashed lines correspond to linear fits with/without 
offset (see insets for goodness-of-fits values). 
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0.5 





Fig. 3. — Probability density functions of the correlation coefficient between 7-ray and X-ray 
fluxes: /a data, fs perfectly correlated case, fc uncorrelated case, and fa analytical solution 
in the uncorrelated case. See text for further details. 
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Fig. 4. — Probability density functions of the correlation coefficient between 7-ray and 
optical fluxes: /a data, fs perfectly correlated case, fc uncorrected case, and analytical 
solution in the uncorrelated case. See text for further details. 



x10- 



H 1.5- 



d 1 



0.5- 



I ' ' ' ' I ' ' ' ' I-' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I 



Crab Nebula 



E 15 



O 5 
U 



I ■ I ■ I I I ■ ■ I n I ■ I ■ I I ■ I ■ I ■ ■ I ■ I ■ I ■ ■ I 



x10'° 



Bkg Rate [counts/min] = 9.9 ± 0.2 (Chl2/NDF = 32.5/30) 
,r1**-««l>i«>«i«i«-iai 



21:50 22:00 22:10 22:20 22:30 22:40 22:50 

Time 



H 1.5 
o 



d 1 



I ' ' ' ' I ' ' ' '- 1 ' ' ' ' I ' ' ' ' I 




Crab Nebula 



I .... I 



il . . . . I 



E 15^ 

I 10 

i 5 

O 



Bkg Rate [counts/minj = 9.4 ± 0.2 (Chi2/NDF = 17.1/21) 



21:50 22:00 22:10 22:20 22:30 

Time 



Fig. 5. — Integrated-fiux LCs of Mrk 501 for the flare nights of June 30 and July 9. Horizontal 
bars represent the 2-minute time bins, and vertical bars denote la statistical uncertainties. 
For comparison, the Crab emission is also shown as a lilac dotted horizontal line. The vertical 
dot-dashed line divides the data into 'stable' (i.e., pre-burst) and 'variable' (i.e., in-burst) 
emission. The horizontal black dashed line represents the average of the 'stable' emission. 
The solid black curve represents the best-fit fiare model (see eq. [2]). The bottom plots show 
the mean background rate during each of the 2-minute bins of the LCs. The insets report 
the mean background rate during the entire night, resulting from a constant fit to the data 
points. The goodness of such fit is also given. 
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Fig. 6. — LC for the night June 30 with a time binning of 4 minutes, and separated in 
different energy bands, from the top to the bottom, 0.15-0.25 TeV, 0.25-0.6 TeV, 0.6-1.2 
TeV, 1.2-10 TeV. The vertical bars denote la statistical uncertainties. For comparison, 
the Crab emission is also shown as a lilac dotted horizontal line. The vertical dot-dashed 
line divides the data into 'stable' (i.e., pre-burst) and 'variable' (i.e., in-burst) emission 
emission. The horizontal black dashed line represents the average of the 'stable' emission. 
The 'variable' (in-burst) of all energy ranges were fit with a constant line. The results of the 
fits are given in the insets. The constant line fit on the energy range 0.25-0.6 TeV was not 
satisfactory (see inset in figure); yet this LC could be fit with the fiare model described by 
equation [2] (see table [3] for the resulting parameters). 
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21:50 22:00 22:10 22:20 22:30 



Fig. 7. — LC for the night July 9 with a time binning of 4 minutes, and separated in different 
energy bands, from the top to the bottom, 0.15-0.25 TeV, 0.25-0.6 TeV, 0.6-1.2 TeV, 1.2- 
10 TeV. The vertical bars denote la statistical uncertainties. For comparison, the Crab 
emission is also shown as a lilac dotted horizontal line. The vertical dot-dashed line divides 
the data into 'stable' (i.e., pre-burst) and 'variable' (i.e., in-burst) emission emission. The 
horizontal black dashed line represents the average of the 'stable' emission. The 'variable' 
(in-burst) of all energy ranges were fit with a flare model described by equation [2], where 
c = d (rise=fall time). All parameters were left free in the fit. All light curves were considered 
simultaneously in the fit (combined fit). The resulting parameters from this combined fit 
are reported in table HI 
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Fig. 8. — Same as in figure [71 but with a common TO (which was also fit) for all LCs. The 
resulting parameters from this combined fit are reported in table [51 
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Fig. 9. — Fractional variability parameter as derived for 5 energy bins. Vertical bars denote 
la uncertainties, horizontal bars indicate the width of each energy bin. The left plot includes 
all data. The right plot includes all but the June 30 and July 9 data. The black horizontal line 
in the right panel results from a constant fit to the data points (see inset for fit parameters). 
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Fig. 10. — Fractional variability parameter as computed for the nights of June 30 and July 
9. Vertical bars denote la uncertainties, horizontal bars indicate the width of each energy 
bin. The black horizontal lines are the result of a constant fit to the data points (see inset 
for fit parameters). 
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Fig. 11. — Correlation between spectral shape and integrated flux above 0.15 TeV. Each 
point denotes a single night of observation. The error bars denote la uncertainties. The 
June 30 and July 9 data were split chronologically into two data sets each, corresponding to 
the pre-flare ('stable') and in-flare ('variable') emission shown in flgureO 




Fig. 12. — Probability density functions of the correlation coefficient between spectral in- 
dex and 7-ray fluxes: /a data, fs perfectly correlated case, fc uncorrelated case, and /d 
analytical solution in the uncorrelated case. See text for further details. 
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Fig. 13. — Measured spectra of Mrk 501 for three different >0.15 TeV flux levels: low (green 
down triangles), medium (red up triangles), and high (black squares). Vertical bars denote 
la uncertainties, horizontal bars denote energy bins. Lines show power-law best fits. See 
table [6] for fit parameters. 




Fig. 14. — Measured spectra for the nights of June 30 and July 9 when Mrk 501 flared. Bars 
as in Fig. [131 Spectral fits are a power-law (red dashed line; see eq. |5]) and a log-parabolic 
function (black solid line; see eq. E]). See table [7| for fit parameters. 
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Fig. 15. — The spectra of Mrk 501 in the nights of June 30 and July 9 corresponding to 
the pre-burst ('stable') and in-burst ('variable') emission (see sect. 13.31 and Fig. E]). Black 
squares/red triangles denote 'stable'/'variable' emission. Bars as in Fig. [121 The insets show 
the log-parabolic fit parameters (see eq. [6]). 





2005 Jun 30 






1 ndf 


4.585 ; 6 




X' 1 ndf 


8.244 ; 8 




Prob 


0.5981 




Prob 


0.41 




pO 


0.1319±0.01948 




pO 


0.1602 + 0.01885 





21:50 22:00 22:10 22:20 22:30 22:40 22:50 

Time 



g 0.7 

8 0.6 
c 

■D 

i5 0-5 

T 

0.4 
0.3 
0.2 
0.1 





2005 Jul 9 






y} 1 ndf 


16.39/5 




x'/ndf 


18.33/5 




Prob 


0.005823 




Prob 


0.002556 




pO 


0.06464 ±0.01983 




pO 


0.1078 + 0.02134 






Fig. 16. — Hardness ratio F{1.2 — 10 TeV) / F(0.25 — 1.2 TeV) vs Time for the nights of June 
30 and July 9. Horizontal bars represent the 4-minute time bins, and vertical bars denote 
la statistical uncertainties. Black squares and red triangles denote pre-burst ('stable') and 
in-burst ('variable') emission respectively (see sect. 13.31 and Fig. [5]). The lines result from a 
constant fit to the data points (see insets for fit parameters). 
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Fig. 17.— Hardness ratio F(1.2 - 10 TeV)/F(0.25 - 1.2 TeV) vs F{> 0.25 TeV) for the 
nights of June 30 and July 9. Horizontal and vertical bars denote la statistical uncertainties. 
Black open squares and red open circles denote pre-burst ('stable') and in-burst ('variable') 
emission respectively (see sect. 13.31 and Fig. [5l). The numbers inside the markers denote the 
position of the points in the LCs. The consecutive (in time) points of the in-burst LC are 
connected by red lines for better clarity. 
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Fig. 18. — Single-night VHE LCs of Mrk 501 obtained with various lACTs during several 
years. Vertical error bars denote la statistical uncertainties. Instruments and corresponding 
mean fluxes are reported for each observational campaign separately. The MAGIC data were 
reprocessed to match the energy ranges covered by previous instruments. 
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Fig. 19. — Spectra of Mrk 501 in the nights of June 30 (left) and July 9 (right) when the 
source was flaring. Bars as in Fig. [131 The spectra have been corrected for EBL absorption 
using (iKneiske et alll2004j )'s 'Low' EBL model. The curves show log-parabolic fits (see eq. 
E]) whose corresponding parameters are reported in the insets. 
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Fig. 20. — Spectral peak location versus flux above 0.15 TeV for the two flaring nights (June 
30 and July 9). The spectra were fitted with eq. [6] (see Fig. [T9l) and the peak location and 
its associated error were calculated using eqs. [7] and [HI The black squares correspond to the 
observed spectra and the blue circles correspond to the spectra after correction for the EBL 
absorption using (iKneiske et al 1 12004 ) 's 'Low' EBL model. 
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Fig. 21. — Overall SED from MrkSOl. The optical data from the KVA Telescope is repre- 
sented with a green full circle; the X-ray data from RXTE / A.'SM. depicted with a black full 
triangle for June 30, red open circle for July 9, and light blue open square for the other 
nights (combined); the VHE data from MAGIC are represented as black full triangles (June 
30), red open circles (July 9), green full circles ('high flux' data-set), dark blue open triangles 
('medium flux' data-set), and light blue open squares ('low flux' data-set). See Sect. 14.21 for 
definitions of high, medium and low flux data-sets. Vertical error bars den ote la statisti- 



cal u ncertainties. The VHE spectra are corrected for EBL extinction using (IKneiske et al 



20041 ) 's 'Low' EBL model. The high e st an d the lowest state were fit with a one zone SSC 



model (described in iTavecchio et al I (120011 )). See table [H] and section [5751 for details of the 
modeling. 



